### Reproduce Appendix Figure 3 (Harris)

### PREP
rm(list = ls())
#setwd(..)
library(foreign)
smoking <- read.dta("harrisxx2.dta")
library(gam)

################################## Figure 3 #############################################

### Relative vs. Attributable risk 
  # Quadrant 2 graph - page 2 of pdf output 
  # Quadrant 3 graph - page 3 of pdf output

smoking.gam1 <- gam(currsmok1 ~ s(rd) + s(rr),
                   subset = control==1&rd>-500,
                   data=smoking, family=gaussian(), weights=weight1)
pdf("harris_attributable.pdf")
  plot(smoking.gam1, xlab="Relative Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
  plot(smoking.gam1, xlab="Attributable Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
dev.off()

### Relative vs. Absolute risk 
  # Quadrant 1 graph - page 2 of pdf output
  # Quadrant 4 graph - page 3 of pdf output 
smoking.gam1a <- gam(currsmok1 ~ s(q2062_) + s(rr),
                   subset = control==1&rd>-500,
                   data=smoking, family=gaussian(), weights=weight1)
pdf("harris_absolute.pdf")
  plot(smoking.gam1a, xlab="Relative Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
  plot(smoking.gam1a, xlab="Absolute Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
dev.off()

